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1 Introduction 

In order to understand the global structure, dynamics, and physical and chemical 
processes occurring in the upper atmospheres, exospheres, and ionospheres of the Earth, the 
other planets, comets and planetary satellites and their interactions with their outer particles 
and fields environs, it is often necessary to address the fundamentally non-equilibrium 
aspects of the physical environment. These are regions where complex chemistry, energetics, 
and electromagnetic field influences are important. Traditional approaches are based largely 
on hydrodynamic or magnetohydrodynamic (MHD) formulations and are very important and 
highly useful. However, these methods often have limitations in rarefied physical regimes 
where the molecular collision rates and ion gyrofrequencies are small and where interactions 
with ionospheres and upper neutral atmospheres are important. 

At the University of Michigan we have an established base of experience and 
expertise in numerical simulations based on particle codes which address these physical 
regimes. The Principal Investigator, Dr. Michael Combi, has over 20 years of experience in 
the development of particle-kinetic and hybrid kinetic/hydrodynamics models and their direct 
use in data analysis. He has also worked in ground-based and space-based remote 
observational work and on spacecraft instrument teams. His research has involved studies of 
cometary atmospheres and ionospheres and their interaction with the solar wind, the neutral 
gas clouds escaping from Jupiter’s moon Io, the interaction of the atmospheres/ionospheres 
of Io and Europa with Jupiter’s corotating magnetosphere, as well as Earth’s ionosphere. 

Grants from the Applied Information Systems Research Program (AISRP) have 
supported the basic construction and first applications of three-dimensional, neutral and ion 
kinetic particle models, which have potential application to various space science problems. 
The purpose of this grant, which resulted from an unsolicited proposal for 1 year of support 
at $58,868, was to expand the original research and to develop the algorithms and operational 
code to enable seamless interaction between kinetic and hydrodynamic descriptions. This 
report describes our progress during the year. The contained in section 2 of this report will 
serve as the basis of a paper describing the method and its application to the cometary coma 
that will be continued under a research and analysis grant that supports various applications 
of theoretical comet models to understanding the inner comae of comets (grant NAG5- 13239 
from the Planetary Atmospheres program). 

2. Hybrid DSMC/CFD Approach 

Simulations of real gas flow often require numerical solutions with a wide variation 
of macroscopic parameters throughout computational domain. In some cases, the gas flow in 
a part of the domain could be in rarefied regime while another part may be in continuous 
regime. Rarefied flow conditions appear in different kinds of circumstances in the numerical 
study of cometary comae (e.g., the Knudsen layer at the surface, on the night side of nuclei, 
at larger distances from the nucleus and for weak comets in general). Even though the kinetic 
approach is generally appropriate physically to study gas flows even in the fluid limit 
(because the fluid equations themselves are approximations of the general kinetic 
formulation), it would make the numerical simulations of the flow very inefficient 
computationally. So, developing of efficient numerical approach that can be used for 
simulation of gas systems under continuum/rarefied condition is an important problem of 
modem numerical gas dynamics. An approach that utilizes DSMC and CFD methodologies 
in different subdomains is described in various papers [Bird, 1998; Ivanov et al., 2000; 
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Ivanov and Markelov , 2000; Hash and Hassan, 1996; Roy et al., 2002; Wang and Boyd, 
2003; Wang et al., 2002]. 


2.1 Dusty-Gas Hydrodynamics 

The appropriate description for the physical state of any one of up to s species, p, in a 
multispecies dilute gas can be described by the Boltzmann equation, which is given as 
(Combi et al. 2004): 

~ ~ ^ s +» 4 * 

-o r )*c r --(f r )*F' •— </,)-2J /(/;/,;-/,/.,> (i) 

where / p s f p (r,c,t) and/ d ■ f d (r,c,t) are the full phase space velocity distribution 

functions for species p and q, cq is the velocity, r is the spatial coordinate, F represents 
external forces, crpq is the magnitude of the relative velocity between particle p and q, & 2 is 
the solid angle, the asterisks (*) indicate post-collision particles, the Iq subscript refers to the 
scattering target particles, and Opq is the total collision cross section between species p and q 
(which can in general be velocity dependent). In the collision integral on the right hand side 
of the equation,/* represents additions of particles, scattered into the region of velocity space 
in question, i.e.. between c and c +dc . whereas /without the asterisk represents scattering 
out of that region. The Boltzmann equation makes neither assumptions nor any restrictions 
as to the form of the distribution functions. The description of the Euler equations for a 
dusty-gas cometary coma have been described in detail in many papers, including some by 
our own group at the University of Michigan ( Combi et al. 1999, Korosmezey and Gombosi 
1989). Various velocity moment expansions of the integrals over the distribution functions 
yield equations of conservation of mass, momentum and energy, which can be written for a 
dusty-gas neutral cometary coma as: 

^ + V-(pu) = ^, 
dt bt 


p— + p(u-V)u+Vp = -F, 
dt 


1 ^ + - 7j( u 'V)p + -^p(V-u) = -(2 g(i +Q ph -Q m 


y - 1 dt y-i y 

respectively, and for any number of dust particle-size populations 


+ v • (p,u,) = i-wv, 

at ot 


Pi + p,(y, • V)V, = F, i = ...,N 

at 

respectively, where p is the gas mass density, u is the gas velocity, p is the gas pressure, v, 
and p, are the velocity and mass density for dust particles of radius a,-. The right hand sides 
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of all the equations contain the various source terms. The term ~ is the gas production 

n n 

source rate, and — is the dust production source rate for particles of radius, a,. F is the gas- 
dt 

dust drag force, which is related to the forces on the individual particle size populations as 

N ~ 

F = -\F ( , where the size dependent force is given by F, = — pC D s i , where 
M 4fl iPw 

U- v 

s,. = , k is the Boltzmann constant, T is the gas temperature and m is the gas mean 

■yjlkTIm 

molecular mass. The accommodation of gas via collisions with hotter dust yields the dust- 
gas heat exchange rate Q g </ which is given as 


o*-—pc r u2ar-w', 

V 
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Here C p is the gas heat capacity at constant pressure and the rest of the coefficients can be 
defined under the assumption of diffusive reflection such that: 
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Finally, pd is the bulk mass density of dust particles of radius a,; T is the gas temperature, T, 
is the dust temperature assumed in equilibrium with solar radiation. Otherwise the other 
intermediate quantities are s„ the relative Mach number between gas and dust, Co the dust- 
gas drag coefficient, Tf ec , the recovery temperature, R,', the heat transfer function and St,', the 
Stanton number. Most of the dust-gas drag physics comes from the formulation by Finson 
and Probstein (1968) with later corrections discussed by Wallis (1982) and Kitamura (1986). 

As described in the paper by Combi et al. (1999), we have successfully applied dusty- 
gas dynamics to a very collisionally thick comet, Hale-Bopp (1995 01). However, for most 
comets the gas becomes too rarefied both for productive comets at large distances from the 
nucleus, as well as for weak comets throughout most of the coma. 
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2.2 Continuum breakdown 

It is well known that there are two major limiting conditions of applicability of the 
computational fluid dynamics (CFD) approach. Fluid description of gas media is valid if 
velocity distribution function of gas molecules is close to Maxwellian and macroscopic 
parameters changes slowly on the scale of mean free path. The physical reasons [Boyd, 1995] 
for breakdown in expanding and compressing flows are different in several ways. In 
expanding flow, continuum breakdown generally occurs when the collision rate becomes so 
low that thermal equilibrium cannot be maintained among the energy modes. In compressed 
flows, continuum breakdown occurs in general due to the very steep flow gradients in the 
shock front, and next to the body surface. 

A continuum breakdown parameter [Wang and Boyd, 2002] is generally used as the 
criterion for switching between methods. In his pioneering work. Bird [1965] first advocated 
a semiempirical parameter for expanding flows. Another empirical parameter based on local 
flow gradients was later developed specifically for hypersonic flows. More recently, a new 
breakdown parameter based on the Chapman-Enskog perturbation expansion of the 
Boltzmann equation has been developed again for expanding flows. The breakdown of 
equilibrium was found to coincide with a certain value of the ratio of the logarithmic time 
derivative of density following the motion of the fluid to the collision frequency of the gas. 
This value had been proposed [Ivanov et al., 1999; Boyd, 1995] as an empirical breakdown 
criterion for use in engineering studies of systems which involve low-density expansion from 
continuum to highly rarefied conditions 
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This form of P may be interpreted as the Lagrangian mean free path divided by the density 
scale length. Breakdown parameter (4) can be evaluated based on other macroscopic 
parameters of the gas system such as mean flow speed V or translational temperature T . In 
this case, the value of the breakdown parameter P is determined [1 1 1] as 


-max(P D ,P r ,P v ) 


( 5 ) 


where P D ,P T ,P V are determined by substituting density number, temperature and velocity 
for p in the expression (4). The breakdown of the continuum assumption usually coincides 
with a degree of translational nonequilibrium for which anisotropic effects arise. Bird [Hash 
and Hassan, 1996] advocates a continuum breakdown parameter (3), with a limit value 
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between 0.01-0.02 

The Knudsen number is the key parameter in determination of the gas flow regime. It 
is commonly assumed [Boyd, 1995] that failure of the continuum approach occurs at 
Knudsen numbers around 0.01. Limiting value of Knudsen number 0.1 Kn is proposed in 
[Hash and Hassan, 1996]. A local Knudsen number may be calculated [Boyd, 1995] 
throughout the flow field based on local flow properties and the length scale of the body. 
This parameter is termed the body-length local (BLL) Knudsen number: 



( 6 ) 


where X is the local mean free path, and D is a fundamental length scale of the problem. The 
limitation of the BLL definition is that a constant length scale, D, is employed. The local 
Knudsen number can better be formulated in terms of length scale determined by flow 
gradients. The gradient-length local (GLL) Knudsen number can be expressed [Boyd, 1995] 
in the following form 


Kn 


GLL ~ 


A 

dQ 

Q 

dl 


( 7 ) 


where A again is the local mean free path, Q is a flow property (density or temperature), and 

/ is some distance between two points in the flow field. This distance should be taken 
approximately along the line of the steepest gradients in the flow properties, streamlines or 
computational grid lines. 

In practical calculations, it is often more convenient to simplify expression (7) by evaluation 
[Wang and Boyd, 2002] dQ/dl as V<2 without projecting it onto a preferential direction. 
Now the GLL Knudsen number for property Q can be written as 

Kn Q - ^|V(2| (8) 

The Knudsen number of this form has a great physical meaning. When its value is much less 
than unity the flow can be regarded as locally slightly perturbed from equilibrium [111] that 
is a fundamental assumption of the fluid approach. A new parameter, max Kn max has also been 
proposed [Wang and Boyd, 2002] that is modified from Boyd's Kn CLL parameter. In the same 
way as it was done in expression (5), the value of the local Knudsen number can be obtained 
based on different macroscopic parameters and define in the form 

Kn mm = max(Kn D ,Kn T ,Kn v ) (9) 

In the rarefied gas dynamics literature it is now generally accepted that continuum 
breakdown is best predicted wherever the value of Kn max exceeds 0.05. 

An entropy generation rate has been proposed [Chen et al., 2003] as a continuum 
breakdown parameter since it has a stronger physical foundation than Knudsen numbers. 
The entropy generation is defined as 


, du , 4/ dT 

gen T dx j T 2 dxj 


( 10 ) 
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Two non-dimensional parameters of entropy generation have been introduced for 
axisymmetric problems: 
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where the former uses local and the latter uses global free stream conditions. The global 
parameter threshold is proposed to be 1, while the local parameter threshold is proposed to be 
5x1 O' 4 

A Grad criterion has been proposed [Tallis and Mallinger, 1997] for determination of 
the Navier-Stokes (NS) solution. For a given known NS solution, moments of the 
Boltzmann equation, K c =(1,£,§-§,|£ 2 |§), need to be evaluated for the corresponding Grad 
Chapman-Enskog distribution function, giving the residual as 

R NS = fK G l^ + y^-Q{f CE )^ , (13) 


where 


e(/./) - / (/r - ff) %|.k ■ H) *» • (W) 

If the NS distribution is a good approximation of the distribution for the general Boltzmann 
equation, the value of the residual Rns will be close to zero, and the NS solution will be valid 
everywhere the residual is small. Similar criteria can be formulated to evaluate the 
discrepancy of distribution function from a Maxwellian which is useful for determination the 
area of validity of Euler equations. 

2.3 The DSMC/CFD Boundary and Interface 

Coupling kinetic and fluid approaches basically involves [Hash and Hassan, 1996] 
determination of the properties at the interface between the DSMC and CFD regions. The 
macroscopic properties are required by the conservation equations for the evaluation of the 
net fluxes and are necessary for the DSMC method for the specification of the distribution 
function used to initialize particles entering the DSMC region. At the end of each hybrid time 
step, particles enter from the fluid region, such that mass, momentum and energy fluxes are 
conserved along the interface. Several methods of maintaining the interface between these 
regions have been proposed. 

Marshak condition 

The most obvious way of coupling DSMC and CFD approaches is based on equating of half- 
fluxes at the interface. The half-fluxes can be defines as the flux of particles which velocities 
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comprising only of the normal velocity component. The DSMC half-fluxes are 


77 j_ num 

DSMC ~~ AM 


my, 

m i v 

— v • v 
2 V 


On the fluid side, the half-fluxes are 

00 00 co 

f cfd = n f f f v Qf du dv d q 


—00 0 -00 


( 15 ) 


( 16 ) 


Q = 


1 

m i v i 



(17) 


The net flux is the sum of the incoming and outgoing half-fluxes. For the Marshak condition, 
the half-fluxes of the mass, momentum and energy of the particles crossing the interface 
during the DSMC step are measured. Theses half-fluxes are then added to the incoming half- 
fluxes from the fluid side and equated to the net flux. When assuming a Maxwellian 
distribution on the fluid side, this gives five equations and five unknowns, providing the 
macroscopic properties of density, velocity and temperature at the interface. 

Extrapolation of flow properties 

This approach of maintaining the interface between DSMC and CFD subdomains 
involves extrapolation of the macroscopic properties of density, velocity and temperature to 
the interface. The accuracy obtained is dependent on the order of extrapolation. The DSMC 
solution is allowed to evolve while accumulating samples for the determination of the 
macroscopic properties [Hash and Hassan, 1996]. After a specified number of time steps, the 
solution is stopped, and the velocity and temperature are evaluated. At the first coupling, 
initial guesses for the velocity and temperature at the interface are taken from the full fluid 
solution. Subsequent initial guesses for the interface values are made using the values from 
the previous coupling. The two solutions are then smoothed across the interface region and 
new values for the interface are obtained. The solution for the fluid region is repeated, and 
the process iterated upon until the interface values converge. Because a good smoothing 
routine is required, a quadratic least-squares fit for the profiles across the interface is 
employed. High-order fits should be avoided because of the oscillations they can produce. 

Extrapolation of the net fluxes 

Instead of extrapolating the macroscopic parameters of the gas flow, net fluxes can be 
extrapolated across the interface. The main difficulty with this is that evaluation of the net 
flux constants must be very accurate in order to ensure that the calculated values conform 


8 


with the profile predicted by the DSMC method. In practice, the scatter in higher moments is 
so great that the method is essentially useless. Example of implementation of the method is 
presented in [Roy et al., 2002.] where a cell-centered finite- volume approach is used for 
discretization of governing equations and second-order reconstructions of the interface fluxes 
is obtained via MUSCL extrapolation. 

Overlapping subdomains 

The approaches considered above are based on separation of the whole computational 
domain into a set of non-overlapping subdomains where an appropriate gas dynamics 
description is applied and where either macroscopic parameters or fluxes of the flow are used 
to couple the subdomains. Another approach, where kinetic and fluid subdomains overlap 
each other has been described in [Boyd et al., 2003]. 
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Figure. 1 Schematic of the interface between continuum and kinetic regions for the case of overlapping 
sub-domains 


Buffer cells and reservoir cells are used in the continuum domain as illustrated in Fig. 
1. These buffer cells and reservoir cells are also treated as particle cells for particle 
movement and particle collisions. Then the interface becomes the internal cell edge for the 
DSMC. The reservoir cells are used to generate particles that can enter the DSMC domain, 
which avoids directly generating particles on the interface. The buffer cells, however, 
improve the quality of the particles that enter the DSMC. As a consequence, Chapman- 
Enskog or Maxwellian distribution functions can be assumed in the fluid region described by 
Navier-Stokes or Euler equations, respectively. Particles are generated in the buffer cells 
according to the distribution function based on the local macroscopic information. In the 
beginning of each time step, particles are generated in the reservoir cells according to the 
distribution function and at the end of the iteration all particles in the reservoir cells are 
removed. Particles in the buffer cells and reservoir cells are selected for collisions and move 
around using general DSMC procedures. The particles leaving particle cells (including buffer 
cells and reservoir cells) are removed. 

Application of the Hybrid DSMC/Euler method to comet 21P/Borrelly 

In September 2001, NASA's Deep Space 1 spacecraft flew by comet 21P/Borrelly 
and obtained only the first set of close-up images and spectra of a cometary nucleus since the 
flyby of comet 1 P/Halley in 1986. Longer exposure images centered on the nucleus revealed 
very interesting, structured sets of dust jets emanating for various source regions on the 
dayside of the nucleus. Our first application for the hybrid DSMC/Euler method is ID- 
spherical gas coma relevant for the nucleus size and gas production rate of comet Borrelly. 
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We first solved for the flow using the Euler solver to get initial conditions for the hybrid 
DSMC/EULER simulation, where the problem is solved in fluid approximation throughout 
whole domain. Even though the value of Knudsen number indicates that gas in some parts of 
the domain is under rarefied regime, this approach allows us efficiently generate initial 
conditions for hybrid simulation which would be “close” to converged solution. It also 
provides an interesting point of comparison to illustrate the weakness of the fluid 
approximation in some regions. To speed up convergence, simulations begin with a coarse 
mesh (i.e., not many computational cells and particles) and subsequently the mesh is refined 
and additional simulation particles added appropriately. 

To couple the DSMC and EULER subdomains, we used interpolation of macroscopic 
parameters on the boundary between subdomains as discussed above. At the start of the 
hybrid simulation, the EULER-solution was used for interpolation. Interpolated values were 
used to generate model particles that enter DSMC region and to calculate fluxes on the 
EULER side, where a standard MUSCL finite volume scheme is used. Model particles on the 
boundary were generated with a Maxwellian distribution function. The main difficulties in 
DSMC/EULER hybrid simulation is due to statistical nature of DSMC method where a 
“large” sample in microscopic parameters is required to get acceptable noise level in the 
macroscopic flow field. In this work, different numbers of iterations were used on EULER 
and DSMC sides. While iterations on the EULER side are required only for convergence of 
the fluid solution, on DSMC side solution is iterated for decreasing of noise level also. 

Table. Physical parameters for Comet Borrelly for the ID spherically symmetrical case. 

Uniform Maxwellian at the surface 
Pure water coma 
Production rate of 10 29 s' 1 
Zero bulk surface velocity 
Surface temperature = 1 80 K 

Model particles that cross boundary of computational domain are deleted. 

Euler domain Kn < 0.03 
DSMC domain Kn >0.03 

A 

Knudsen number in the form Kn n = — 

2 Q 

Figure 2 shows our first results using the new hybrid model compared with Euler 
hydrodynamics alone. The most obvious differences are (as should be expected) in the 
expanding coma. In the hybrid model the isotropic Maxwellian distribution is clearly not 
present, as evidenced by the divergence of the radial and surface-normal components of the 
temperature. As in the standard rocket nozzle problem, the radial temperature freezes out at 
a non-zero value. The additional effect of this is that the asymptotic radial expansion 
velocity is smaller in the hybrid model than in the Euler calculation. This is caused by the 
fact that the radial temperature (i.e., the radial dispersion in velocities) is not transferred to 
directed radial motion. 

Another interesting result is subtler. In comparing a number of curves (but 
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Comparison of Pure Euler Hydrodynamics with Hybrid Euler/DSMC. On the left is a model calculation 
using a standard 5 -moment Euler hydrodynamics scheme for the expansion of the cometary atmosphere from 
the surface of the nucleus (5 km) to about 100 km (i.e., 10 7 meters). The surface Knudsen layer is assumed to 
be a sonic transition (a shock) which is solved using the Riemann conditions. On the right the hybrid scheme is 
used. The thin Knudsen region near the nucleus (between the y-axes and the tall rectangular boxes is solved 
using DSMC. The rectangular box is the fluid region (Kn < 0.03), and to the right of the box is a second kinetic 
transition region. Clearly the temperature plot (middle panel right side) shows the separation of radial (top), 
isotropic (middle) and tangential (bottom) components owing to fluid breakdown. Density is in units of m* 2 , 
temperature in Kelvins and velocity inms' 1 . Distance (x/R is normalized to the nucleus radius R=5 km. 
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especially the density and temperature profdes) at the outer boundary of the Euler region in the 
hybrid calculation (this would be at the right edge of the rectangular box on the plots on the left 
hand side of the figure), it appears that although the macroscopic value is continuous across the 
Euler-DSMC boundary, the slope is slightly discontinuous. This is either very interesting (or 
perhaps troubling), but in any worthy of further study. Although an error in either the DSMC or 
Euler solver (or in their interface) is always a possibility, it is worth stressing that each code (the 
MUSCL finite element solver for the Euler equations and the DSMC method) has been 
rigorously test against a number of classical problems where this type of behavior did not occur. 
In testing of our first DSMC code [Combi, 1996] we did not run enough simulation particles to 
reduce the noise enough to see such a small and subtle change. An alternative possibility is that 
for this type of expansion problem (i.e., free expansion into a vacuum) even the very 
conservative boundary at Kn = 0.03 is not small enough and that the continuous expansion 
means we have to move the boundary to even smaller values of Kn. 

We will continue to study this problem along these lines with funding from a comet data 
analysis grant which supports the application of models developed here. As mentioned already 
this grant will support completion of a model and comet application paper which will contain 
much of the text that is part of this report. A number of future studies will be to model a much 
smaller Knudsen number regime ( 0.01 or less) in which one would think the DSMC and Euler 
would provide exactly the same results and then see if it is truly the extreme nature of this 
physical problem which is causing the breakdown. 

Finally, one of the important goals of the AISRP is to provide new computational tools 
for use in future NASA programs. An interdisciplinary international modeling effort was begun 
this year, organized by the International Space Science Institute (ISS1) in Bern, Switzerland, in 
association with the Rosetta project (ESTEC) and the JPL Rosetta project (the support and 
organizing infrastructure for US investigators on the Rosetta project). As part of this effort, this 
PI is the lead investigator in the area of dusty-gas coma modeling (with T. Gombosi, U. 
Michigan and N. Thomas, U. Bern). The larger effort involves several other teams of 
investigators performing (1) MHD modeling of the comet/solar-wind interaction (T. Gombosi, 
K. Hansen, and M. Combi, U Michigan), (2) nucleus and surface-boundary modeling (B. 
Davidsson, ESTEC and H. Rickman), (3) electron energetics (T. Cravens U. Kansas), (4) plasma 
hybrid modeling (U. Motschmann and T. Bagdonat, MPI Braunschweig), (5) coma chemistry (K. 
Altwegg, U. Bern, and the entire U. Michigan group). Rosetta Project Scientist G. Schwehm and 
US Rosetta Project (JPL) Scientist C. Alexander are also participants as are a number of other 
Rosetta scientists. The hybrid Euler/DSMC model will likely be an integral part of this modeling 
effort for the next 12 years or more. 
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